微生物组学数据中的污染物的识别与去除 | 微生物专题

在测序过程中扩增得到的有些DNA序列实际并没有真正出现在样本中，污染可能包括方方面面，例如测序反应中使用的试剂、实验操作问题等，可能干扰下游分析，错误的扩大了样本的多样性，影响了数据的准确性，掩盖了正常的数据。虽然适当的实验室操作可以减少污染，但不能消除污染。

这里我们介绍Decontam (https://github.com/benjjneb/decontam)。它是一个开源R包，通过统计分类过程来识别污染物和去除污染的DNA序列，提高测序结果的质量。

原理:

1.   污染物在低浓度样品中出现频率较高；

2.   污染物出现在阴性对照品中出现频率较高

数据准备：

1.   特征表：每个样本(行)中序列特征(列)的相对丰度表；这些“序列特征”可以是各种各样的特征类型中的任何一种，包括扩增子序列变体(ASVs)、操作分类单元(OTUs)、物种分类等。

2.   用于获得污染参考的数据，以下二选一

（1）DNA quantitation data：包含每个样品中DNA浓度的DNA定量数据文本。

（2）negative control：特征标中存在没有添加任何生物样本的空白上进行测序的样本数据。值得注意的是阴性控制也应通过PCR步骤进行，因为工作流程中的每一步都有可能引入新的污染物。

具体步骤：

1.   Setting up

先导入R包及表格

library(phyloseq)library(ggplot2)library(decontam)

ps <- readRDS(system.file("extdata", "MUClite.rds", package="decontam"))

Ps

## phyloseq-class experiment-level object

## otu_table() OTU Table: [ 1951 taxa and 569 samples ]

## sample_data() Sample Data: [ 569 samples by 6 sample variables ]

## tax_table() Taxonomy Table: [ 1951 taxa by 6 taxonomic ranks ]

quant_reading：通过荧光强度测量每个样本中的DNA浓度，

Sample_or_Control：阴性对照情况说明。

head(sample_data(ps))

##                    X.SampleID quant_reading        Sample_or_Control

## P1101C01701R00 P1101C01701R00          2829      True Sample

## P1101C01702R00 P1101C01702R00          5808      True Sample

## P1101C01703R00 P1101C01703R00          5052      True Sample

## P1101C08701R00 P1101C08701R00          1227      True Sample

## P1101C08702R00 P1101C08702R00          3872      True Sample

## P1101C08703R00 P1101C08703R00          5872      True Sample

2.Identify Contaminants - Frequency

2.1 方法一 “frequency”method

该方法利用每个序列特征的频率分布作为输入DNA浓度的函数来识别污染物。

在我们的phyloseq对象中，”quant_reading”一列的记录DNA浓度信息

head(sample_data(ps))

##                    X.SampleID quant_reading        Sample_or_Control

## P1101C01701R00 P1101C01701R00          2829      True Sample

## P1101C01702R00 P1101C01702R00          5808      True Sample

## P1101C01703R00 P1101C01703R00          5052      True Sample

## P1101C08701R00 P1101C08701R00          1227      True Sample

## P1101C08702R00 P1101C08702R00          3872      True Sample

## P1101C08703R00 P1101C08703R00          5872      True Sample

使用isContaminant函数进行污染物的识别，method填写为frequency。

p列表示属于污染物分类的概率，

contaminantl 列的 TRUE/FALSE表示序列特征是污染物的统计值是否超过了用户可设置的阈值，如果为True认为该特征序列为污染物序列（默认阈值为0.1）。

contamdf.freq <- isContaminant(ps, method="frequency", conc="quant_reading")

head(contamdf.freq)

## freq prev p.freq p.prev p contaminant

## Seq1 0.323002694 549 1.000000e+00 NA 1.000000e+00 FALSE

## Seq2 0.098667396 538 1.000000e+00 NA 1.000000e+00 FALSE

## Seq3 0.003551358 160 1.135975e-18 NA 1.135975e-18 TRUE

## Seq4 0.067588419 519 9.999998e-01 NA 9.999998e-01 FALSE

## Seq5 0.045174743 354 1.000000e+00 NA 1.000000e+00 FALSE

## Seq6 0.040417101 538 1.000000e+00 NA 1.000000e+00 FALSE

1901条序列中有58条被认为是污染序列

table(contamdf.freq$contaminant)

## ## FALSE TRUE

## 1893 58

下图给我们可视化呈现了非污染序列Seq1与污染序列Seq3的区别.

在这个图中，横坐标表示DNA浓度，纵坐标表示特征序列的频率。黑色虚线显示了一个非污染序列特征的模型，该特征的频率预计与输入的DNA浓度无关。红线表示的是一个污染物序列特征的模型，该特征的频率预计与输入DNA浓度成反比，因为在DNA总量很少的样本中，受污染的DNA将占总DNA的更大比例。很明显，Seq3非常适合红色污染物模型，而Seq1不适合。

plot_frequency(ps, taxa_names(ps)[c(1,3)], conc="quant_reading") + xlab("DNA Concentration (PicoGreen fluorescent intensity)")

把它们从phyloseq对象中移除

ps.noncontam <- prune_taxa(!contamdf.freq$contaminant, ps)

ps.noncontam

## phyloseq-class experiment-level object

## otu_table() OTU Table: [ 1893 taxa and 569 samples ]

## sample_data() Sample Data: [ 569 samples by 6 sample variables ]

## tax_table() Taxonomy Table: [ 1893 taxa by 6 taxonomic ranks ]

2.2  方法二 “prevalence” method

在这种方法中，将真阳性样本中各序列特征与阴性对照中的比率进行比较，以确定污染物。

在我们的phyloseq对象中，“Sample_or_Control”是包含阴性样本的信息，Control Sample表示阴性样本，True Sample表示阳性样本。

head(sample_data(ps))

##                    X.SampleID quant_reading        Sample_or_Control

## P1101C01701R00 P1101C01701R00          2829      True Sample

## P1101C01702R00 P1101C01702R00          5808      True Sample

## P1101C01703R00 P1101C01703R00          5052      True Sample

## P1101C08701R00 P1101C08701R00          1227      True Sample

sample_data(ps)$is.neg <- sample_data(ps)$Sample_or_Control == "Control Sample"

使用isContaminant函数进行污染物的识别，method填写为prevalence。

contamdf.prev <- isContaminant(ps, method="prevalence", neg="is.neg")

这种方法识别出的污染序列数量为124，比第一种的方法识别出的污染物数量为58要多，但漏掉了污染物Seq3，主要是因为Seq3几乎存在于阴性和阳性等所有样本中。

table(contamdf.prev$contaminant)

## FALSE TRUE ## 1827 124

与第一种方法一样，污染物的识别默认阈值是0.1的概率，我们可以选择尝试将threshold改成0.5。

contamdf.prev05 <- isContaminant(ps, method="prevalence", neg="is.neg", threshold=0.5)

此时，识别出的污染序列数量增加至142个

table(contamdf.prev05$contaminant)

## ## FALSE TRUE ## 1809 142

可视化展示下序列在阴性对照和阳性样本中被观察到的次数，看看这个分类阈值是否可行。
ps.pa                  <- transform_sample_counts(ps, function(abund) 1*(abund>0))

ps.pa.neg          <- prune_samples(sample_data(ps.pa)$Sample_or_Control == "Control Sample", ps.pa)

ps.pa.pos    <- prune_samples(sample_data(ps.pa)$Sample_or_Control == "True Sample", ps.pa)

df.pa          <-   data.frame(pa.pos=taxa_sums(ps.pa.pos), pa.neg=taxa_sums(ps.pa.neg), contaminant=contamdf.prev$contaminant)

ggplot(data=df.pa, aes(x=pa.neg, y=pa.pos, color=contaminant)) + geom_point() + xlab("Prevalence (Negative Controls)") + ylab("Prevalence (True Samples)")

样本分成了两个分支，一个分支主要出现在阳性样本中，另一个分支主要出现在阴性对照中。看来默认概率阈值在识别阴性对照中的污染物方面也已经做得很好了。

参考文献：

Simple statistical identification and removal of contaminant sequences in marker-gene and metagenomics data Davis et al. Microbiome (2018) 6:226

http://www.bioconductor.org/packages/release/bioc/vignettes/decontam/inst/doc/decontam_intro.html

相关阅读